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Mathematical Modeling describes a process and an object by use of the math- 
ematical language. A process or an object is presented in a “pure form” in 
Mathematical Modeling when external perturbations disturbing the study are 
absent. Computer simulation is a natural continuation of the Mathematical 
Modeling. Computer simulation can be considered as a computer experiment 
which corresponds to an experiment in the real world. Such a treatment is 
rather related to numerical simulations. Symbolic simulations yield more than 
just an experiment. They can be considered as a transformation of a mathe- 
matical model by computer, since symbolic simulations keep parameters of the 
model in symbolic form that corresponds to a set of actual experiments. One 
can obtain numerical results as in actual experiments only after substitutions 
of the symbolic parameters with the numerical data. Therefore, symbolic sim- 
ulations complete the mathematical model and embrace actual experiments. 

Mathematical Modeling of stochastic processes is based on the probabil- 
ity theory, in particular, that leads to using of random walks, Monte Carlo 
methods and the standard statistics tools. 

Symbolic simulations are usually realized in the form of solution to equa- 
tions in one unknown, to a system of linear algebraic equations, both ordinary 
and partial differential equations (ODE and PDE). Discrete methods such as 
Finite Element Methods or Finite—Difference Methods for finding approxi- 
mate solutions of ODE and PDE are usually used in the form of numerical 
simulations. 
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1.1 How to develop a mathematical model 
1.1.1 A simple mathematical model 


Consider a simple example of the free falling object from a height h using 
Mathematical Modeling. How to describe the trajectory of the object? The 
first questions which should be posed are the questions “where” and “when”, 
i.e., we have to describe the space where the free fall is taking place and the 
time when it happens. Let gravity be the only force acting upon the falling 
object. In this case, only one axis is needed to describe the space. Let it be the 
axis OY shown in Fig.1.1. It is necessary to fix a unit segment on the axis which 
shows the direction and the line unit. In Fig.1.1, the axis OY is chosen to point 
in the downward direction. Meters can be chosen as the length unit on OY. 
However, dimensionless units are frequently used in Mathematical Modeling, 
which enables the result to be obtained in the easiest way. Dimensionless 
results are transformed into dimension form at the last stage. We will discuss 
this question later in Sec.1.4. Now, we are going to fix a dimensionless line 
unit. The time unit t is also fixed as a dimensionless unit. It is assumed that 
the object begins to fall at the initial time t = 0. 

After answering the questions “where” and 
“when” we should analyse conditions of the 
fall. If the object is a stone, the air resistance 
can be ignored. However, it may not be done 
for a leaf. Let us study a stone now. Then one 
can consider the stone as a material point, i.e., 


y@ as a point on the OY axis for each time t. 
Next, one should use a physical law of the 
falling material point. It follows from the grav- 
itational law that an object near the surface of 
the Earth falls with the constant acceleration 


m 


g (equal to 9.815 in the SI units). Projectile 


" motion is a motion in which a material point 
moves along a trajectory under the action of 
gravity only. Such a material point is called 

Jy projectile . 

FIGURE 1.1: Falling of ma- Mathematics, “hard artillery”, arises at the 


following step to formally describe reality. We 
are interested in representing the motion of the 
point in time as a function y(t) which shows the coordinate of the material 
point in time t. The function y(t) is called the trajectory in mechanics. The 
velocity of the point is the derivative with respect to time v(t) = y’(t), accel- 
eration is the second derivative a(t) = y(t). Now, one can write the law in 


terial object (point) 
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the mathematical form 


“point falls with constant acceleration g” a(t) =g y(t) =g (1.1) 


We have got the ordinary differential equation (1.1). Its general solution has 


the form 
gt? 
y(t) = “+ Cit + Co, (1.2) 


where Cı i C2 are arbitrary constants. The solution (1.2) of equation (1.1) 
is obtained by double integration of (1.1), where Cı and C2 are constants of 
integrations. The trajectory of the material point has to be uniquely deter- 
mined and could not depend on the arbitrary constants C1 and C2. How to fix 
them? First, the law (1.1) does not take into account the information about 
when and how the stone is thrown. In our case, the material point is located 
at the point y = 0 of the axis OY at the time t = 0. Hence, 


y(0) = 0. (1.3) 


Moreover, the object can be thrown with any initial velocity vo. If one just 
drops the stone, the initial velocity is equal to zero. If one applies a force and 
throws the stone down or up the initial velocity vo does not vanish. Hence, we 
have the second condition 

y'(0) = vo. (1.4) 


The conditions (1.3)—(1.4) are called initial in the theory of ordinary differ- 

ential equations and the problem (1.3)—(1.4) is called Cauchy’s problem. This 

problem is properly posed and has a unique solution. It is necessary to substi- 

tute C1 = vo i C2 = 0 in the general formula (1.2) in order to fulfil the initial 

conditions (1.3)-(1.4). As a result we obtain the required particular solution 

2 

y(t) = F + vot. (1.5) 

The mathematical problem has been solved. The height h has to be taken 

into account to find the time T when the projectile reaches the earth. The 

following equation with respect to T is obtained by substitution of t = T and 
y = h into (1.5) 

gT? 


Out of two solutions 


T= —vo + \/2gh + vg T = —vo — y 2gh + vg 
g J 


we must take the first one, since the second solution is always negative which 
contradicts the physical treatment of the problem (T > 0). 
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Remark 1.1. It is a common practice to precisely describe “operating condi- 
tions” in pure mathematics. For instance, in order to work with a function 
mathematicians must define a domain of the function. If a derivative is used, 
a mathematician must say “let a function y(t) be differentiable in a set A”. If 
a mathematician writes “, before this he/she must write m 4 0. However, in 
applied mathematics the mass m is always non-negative. Thus, < for m = 0 
has no any sense. For instance, an applied mathematician writes y’(t) assum- 
ing that this derivative exists. So, usually in applied mathematics if something 
is written, it is tacitly assumed that it exists. 


1.1.2 Use of a computer 


The above example has been solved “by hand” without any computer, 
since the calculations are simple. It is in accordance with the principle: 


Principle of hand calculations. If you can calculate something without 
any computer, try to do it “by hand” If you cannot or you are not sure in 
your calculations, try to use a computer. 


For instance, let somebody not be sure in equations (1.5)—(1.6) based on 
derivatives and integrals. Let us solve the problem by use of the operator 
DSolve in Mathematica! 


Inft;= DSolve[{x''[t] == g, x[0] = 0, x' [0] = vO}, 
x[t], t] 


Out[1}= { {x(t s (gt? +2tvo)}} 


One can copy the resulting function and further paste it into the definition 
of the function y(t). However, it is easier to extract y(t) from the expression 
Out[1] as follows 


infzj:= y[t_] = %[l, 1, 2] 


1 2 \ 
ouia > (gt? +2+t v0) 


Here, % means that the latest result is taken, i.e., Out[1], the triple 1,1, 2 
denotes the coordinates, so 1,1 yields the expression? Rule[x[t], + (gt? + 
2tv0)]. Then, 2 means the second element out of this expression. 

Solve the equation (1.6) and take the proper solution: 


lIn Mathematica, vo means the value of the function v at the point 0. We write the 
symbol v0 instead of vg in Mathematica’s code in order to use v0 as a separate independent 
variable. 

?You can investigate full form the expressions with the FullForm operator. 
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Inf3]:= Solve [h == y[t], t] 


SS -v0+./2gh+v0? }} 


Out[3]= {{e F ' {t 3 


inaj:= T =%[2, 1, 2] 


-v0 +4/2 gh + v0? 


g 


Out[4]= 


Fix the values of the parameters g, h, vo and prepare a graph to demonstrate 
the location of the material point in time. It can be done by use of the operator 
Plot as Plot[h-y[t], {t,0,T}]. To get a better visualisation one can 
do it in a more complicated way: 


In[5]:= {g = 9.81, h=0, vO = -2}; 


Inf6|:= Show[Plot[h-y[t], {t, 0, T}, AxesLabel > {"t", "z(t)"}, 
LabelStyle > 16, PlotRange > {{-0.02, 0.45}, {-0.03, 0.25}}, 
PlotStyle -> Directive [Darker [Gray, 0.1], Thick], 
AspectRatio > Automatic, TicksStyle > 12], 
Graphics [ {Text [Style["T", FontSize > 15, Bold], {1.03T, -0.015}], 
Text [Style["0", FontSize > 15, Bold], {0.03T, -0.015}], 
{PointSize [Medium], Point /@ {{T, 0}, {0, 0}}}}], 
ImageSize > 450] 
z(t) 
0.257 


MATLAB Example Box 1.1 


When computations are performed with the MATLAB package, one 
can use Command Window to enter commands at the command line 
(indicated by the prompt »). Let us define variable g 
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>> g = 9.80665; 


Now, we transform the considered equation (1.1) into a system of ODEs. 
Using yi, = y*-)) (k = 1,2,3) we obtain 


with the corresponding initial values 
yı (0) =h, 
y2(0) = vo. 


Let us define the initial values and place them into a single array 


2> i = =} 
=> w = =a 
>> init_vals = [h, v0] 


Define the function representing right sides of the above system 


F(t, y1, y2) = (yo(t), 9) (1.7) 


using anonymous function as follows 

>> £= @(t, y) ly(2)7 gl; 
We are going to solve the equation numerically, hence we have to provide 
the range of the integration 


>= č imtarval = otis 


Now we are ready to use the built-in solver function ode45 


>> ll, “i oden S (E E limte rva l amni vals); 


Here, the solver returns two arrays: the array T contains time points 
corresponding to rows of the solution array Y. One can combine both 
arrays and plot the solution 


=> plore =s 2) 


The other way of computing in MATLAB, adapted in this book, is 
running scripts written in the MATLAB Editor (files with code organized 
by functions saved with .m extension). Below the source code of file 
script01.m is presented. Note that the very first function, which is 
called the primary function, should have the same name as the name of 
the file. 


% primary function 


function script] () 


@ initial values 
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wo = S26 
in = =2ip 
toincerveal = [0 1)? 


ution of the ODE 


S J 
[T, Y] = clots (6f; t_interval, [h v0]}; 


S SO 


6 plot 

PE IL): @ height of the body 
PLENI H) 

xlabel ( ) 

ylabel ( ) 

grid on 


ylim( [0 max {1 .1*H))) 


@ subfunction 
function dy = f(t, y) 
g = 9.80665; 


dy = zeros(2,1); % a column vector 

e yk denotes (k-1)st derivative of y 
dy(1) = y(2); ® yl’ (t) = y2(t) 

dy (2) = g; yA EES 


The script can be run simply by pressing the Run button from the Editor 
as well as by typing its name into the command window 


== scriptol 


3.5 F T 7 


2.57 ~~ J 


height 


1.57 h3 q 


0.5 + \ 4 


0 L L L L 1 1 L ji L 
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 


time 


The graph of the function z(t) = h— y(t) is drawn with all the numerically 
given parameters. After such a success one can think that an ordinary dif- 
ferential equation can be solved by means of computer without any course of 
calculus. It is enough just to write DSolve or something else in another pack- 
age, to press a button and everything is displayed. But the question is, which 
buttons should be pressed? A computer is a continuation of our knowledge. 
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In order to perform simple computations and to solve simple problems a user 
need not know calculus. But a step aside from a standard way could lead to 
trouble. A poorly prepared user could rather repeat manipulations but could 
not solve a new simple problem. It is possible to give many examples when 
incorrect or superficial use of computer by engineers yields incorrect results. 
Below, we present such an example. 

Calculate 10 partial sums of the series beginning with the partial sum 1160 


for series 
loc) 
k=1 


We have? A E EE IAE E A E E E E EE 


(1.8) 


w| = 


k=1 


In[2]:= Table[S,, {n, 1660, 1670}] 


Out[2]}= {7.99209, 7.99269, 7.99329, 7.99389, 7.9945, 
729951, 7-9957; 7.9963, 7.9969, 7:9975; 7-998097 


Direct observations of the results could yield the conclusion that this series 
converges to 8. However, this series (1.8), which is known as the harmonic 
series diverges*. In order to avoid such mistakes one should use the following: 


Principle of the stupid computer. Do not trust the computer. If you 
do not know what are you computing, use the weaker principle of hand calcu- 
lations. If it is possible to check the computer result “by hand”, do it. 


However, the customer is not always right. A computer gives the correct 
answer to a question, but the user does not understand it and blames the 
computer on the basis of the above principle. Here, the following principle has 
to be taken into account: 


Presumption principle of innocence of the computer. The answer 
you get depends on the question you ask. 


Presumption principle of innocence of the computer completes the prin- 
ciple of the stupid computer. The computer answers the question and the 
computer does not care that the question is stupid. Its task is to answer the 
question in general, without comment, as in the army. It is worth mentioning 
that Mathematica sometimes warns us. Below, Mathematica warns about 
the divergent series 


3The symbol n_ means for all n (Vn) in Mathematica’ s codes; N_ or X_ must be used 
for any argument in definitions of functions. 
ĉit is possible to say in wider sense that it converges to +00 
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<21 
In[3]:= > =i 
k=1 k 


: Sum does not converge. >> 
and division by zero 


ma t1 
n= 


1 
: Infinite expression — encountered. >> 
0 


MATLAB Example Box 1.2 


The following script plots partial sums of another divergent series 


K 1 
BTO (1.9) 


function scripto0z () 

cay MOG Ipdexes fA Se ee -r 2000) 

% corresponding sequence and partial sums 
sequence = 1.0 ./ (k .* log(k)); 
partial_sums = cumsum(sequence) ; 

& plot 


lle (iz, parn iia ESUS, ) 
grid on 


2.5 r 


05 fi 1 i fi 1 i 
0 200 400 600 800 1000 1200 1400 1600 1800 2000 
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As we can see, such visual observation may also suggest convergence of 
the series. Note that the operators ./ and . * are examples of generalized 
arithmetic operations. Here, division and multiplications are performed 
element by element resulting in a new array. Here, the cumsum operator 
computes the cumulative sum (a vector of partial sums of a given vector). 


1.1.3 Development of mathematical models 


In Sec.1.1.1, steps of the development of a mathematical model are demon- 
strated in a simple example. These steps can be applied to general models. 


Steps taken to create a mathematical model: 


i) to introduce spatial variables (description of geometry) and time; 

ii) to think about dimensional units; to introduce the units perhaps during 
solution of the mathematical problem; 

iii) to introduce assumptions and conditions, first, as simple as possible; 

iv) to formulate the law (physical, economical, biological, empirical etc.); 

v) to develop a mathematical description, first, as simple as possible; 

vi) to try to solve the mathematical problem “by hand”; if that does not 
work, to try a numerical method; to compare the results if different methods 
are applied; 

vii) verification of the model; if the results are suspect to get back to the 
previous steps; 


The following comment on the above steps is important. 


Principle of the simplest model. At the beginning do it in the simplest 
way, even it does not yield acceptable results.° 


In the step vi), it is worth doing most of the work “by hand” in order 
to simplify computations. We can say even more. For instance, instead of the 
computation of the series S = ae + it is worth looking up in the handbook 
(for example, implemented in Mathematica) and rewriting the result S = T 

Thus, it is not necessary to perform calculations “by hand”®. For instance, 
the best method to calculate an integral is to find the answer in the hand- 
book’. Hence, the most important is to know where to find the answer. This 
requires good mathematical education. To know where to look is also a type 


of knowledge. In this case, another type of memory is required, not the one 


5¢.f. Albert Einstein: “Everything should be made as simple as possible, but not simpler.” 

Sif only it is useful as an exercise 

"Information as a part of knowledge becomes very important in the modern world. An 
apt example is Google Search used by students and researchers more often than printed 
textbooks 
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that is frequently required to pass exams. In particular, development of math- 
ematical models means the proper selection of existing models scattered in 
textbooks, internet etc. 

During development of the model one should mind that the model has 
to correspond to the real phenomena both qualitatively and quantitatively. 
This contradicts the principle of the simplest model. Anyway, we recommend 
following this principle gradually tending to be in correspondence with reality. 

After verification of the simple model, this model can be improved by 
advanced study. For instance, if instead of the stone (see Sec.1.1.1) a leaf is 
falling, it is better first to discuss the falling of the leaf in vacuum. Then, it is 
possible to complete the model taking into account air resistance. Following 
this approach we pass from the stone to the leaf below. The equation of the 
force equilibrium is similar to (1.1): 


my"(t) =m qg, (1.10) 


where the exterior force ma(t) = my” (t) presented in the left part in accor- 
dance with Newton’s law. The gravitational force mg is in the right part. The 
corrected model for the leaf has the form 


my" (t)=mg—k v(t), (1.11) 


where F = —kv(t) with k > 0 is the air resistance. It is assumed here that the 
air resistance is proportional to the velocity of the fall v(t). It is an empirical 
law, sometimes it has the form F = —kv?(t) or may be more complicated. 
Then equation (1.11) becomes 


V+ yO =g. (1.12) 


The differential equation (1.12) with the initial conditions (1.3)-(1.4) has a 


k 
h= DSolve[{y''[t] = g-—y'[t], y[0] = 0, y'[0] = vo}, yit], t] 
m 


Et kt kt kt 
e = m|gn-e7 gm+em gkt-kvọo+em kvo 
Out[1]= {{y[t: > 5 
k 
Inf2):= %[1, 1, 2] // Expand 
ake kt 
gm? em gm? gmt mvo e m MVọ 
Out[2]= - + + + E 
k? k? k k k 


The trajectory y(t) has been constructed. Further investigation and compu- 
tations can be performed by the scheme presented in the previous Secs 1.1.1- 
112. 


8Zvi Artstein [9]: “Nature is a very good approximation of Mathematics”. 
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MATLAB Example Box 1.3 


One can simulate the considered model numerically as well and ob- 
serve how the acceleration decreases due to the air resistance. Consider 
the following MATLAB script. 


function script03() 


% parameters 
g = 9.80665; 
v0 = 0; 

m = ip 

k = 1; 

h = 100; 


t_interval = [0 10]; 
lution of the equation 


7 
f= el (ie, y) f(t, Yr vO, m, k, g); 
x] = odë45(f; t_interval; [0 v0]); 


i = joy = 2S, ip 
Vv 


I 

l 
x< 
N 


subp roti P2 i) 5 
PLOE(T HI 
xlabel ( ) 
ylabel ( ) 
grid on 


subplot (27.2); 3 aeie subplot 
plot (T, V) 

xlabel ( ) 

ylabel ( ) 

grid on 


UNCclen chy = ee, vr WO, Wy, kp i) 
dy = zeros(2,1); 
dy(1) = y(2); 
dy (2) = g - k/m x y(2); 
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-6t | 
\ 
40 J \ 
A 
L | \ 
30 B r 
20 J -9 + \ 
a 
10 -10 + 
(0) 5 10 0 5 10 
time time 


The final solution should not be given by exact formulae. Consider the 
example of the cooling coffee described in the excellent book by Amelkin [5]. 


Example 1.1. Anatol and Vladimir ordered two cups of coffee with cream in 
a bar. Two cups of hot coffee with the same temperature are simultaneously 
prepared. Anatol and Vladimir take the following steps. Anatol adds cream 
to coffee and covers the cup by a napkin. Vladimir covers the cup by a napkin 
and adds cream to his coffee 10 min later when they simultaneously begin to 
drink their coffee. Who drinks hotter coffee? 


In order to solve this problem in accordance with the book [5] one can 
develop a mathematical model of the cooling coffee based on the theory of 
heat conduction, solve the stated problems and calculate the temperature of 
Anatol’s and Vladimir’s coffees 10 min later. 

However, in this case we are interested only in the question, in which cup 
the coffee is hotter. We are not interested in the values of the temperature. 
The stated problem can be formulated in another way. When does coffee lose 
less heat in 10 min? Does it happen when Anatol decreases its temperature 
by adding cream at the time t = 0 min or when coffee is cooled at the time 
t = 10 min? In order to answer this question we need only one law of the 
heat conduction: dissipation (loss) of heat is proportional to the difference of 
temperatures of the object (coffee) and environment. Hence, hotter coffee loses 
more heat. Therefore, Anatol drinks hotter coffee than Vladimir. 

Similar questions arise in other regimes of heating and cooling, for instance 
in buildings. In order to reach a given temperature at a fixed time with minimal 
dissipation one has to heat the object at that fixed time, not earlier. 

The above problem is based on the estimations of energy. We do not take 
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into account the daily change of prices of energy’. 


Principle of the energy estimations. Estimate energy before you make 
a deal. 


This principle helps us to answer the following questions. Is the force of the 
artillery gun sufficient to overcome the gravitational force to go into space? 
(the answer is negative) Can the Moon cause an earthquake? “Does the flap of 
a butterfly’s wings in Brazil set off a tornado in Texas?”!° Edward N. Lorenz 
asked this question and illustrated it by an example of dynamical system 
which describes the butterfly effect when a small change in one state of a 
deterministic system can result in large differences in a later state. Such a 
question is closely related to stability discussed in Sec.1.3. 


1.2 Types of models 


Types of the mathematical models can be the following: 


a) deterministic and stochastic (random); 
b) discrete and continuous; 
c) linear and non-linear. 


a) An example of the deterministic model is given in Sec.1.1.1, see for 
instance equation (1.12). Let one check the model by means of experiments and 
see that a light wind impacts the result. But how to introduce this wind into 
the model? How to estimate the random influence of the wind? Equation (1.12) 
was useful but should be revised. Let us consider the question from another 
point of view. Introduce the wind into equation (1.12) by force, namely, add 
a random value €(t): 


my" (t) = mg — ky(t) + m€(t). (1.13) 


If the vertical component of the wind force m&(t) is known, the differential 
equation (1.13) can be solved. However, €(t) is just a symbol which does not 
solve the problem and hides it in a new symbol. This term (t) will be useful 
if it is described by the rules used in the probability theory. Below, we outline 
a simple scheme based on the standard course of probability. 

Following the description of the normal distribution we can describe 


°The cost of electric energy at night is usually lower than during day. Moreover, the cost 
at night can be negative. Hence, the economic problem to minimize costs complicates the 
energy problem. 

lWnttp://eaps4.mit.edu/research/Lorenz/Butterfly\_1972.pdf 
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the random function E(t) by two constants, the mathematical expectation 
a = E[€(t)] and the standard deviation o = \/E[E?(t)] — E?[E(t)]. In the 
considered case, a expresses the expected averaged wind force acting on the 
body at time t. Let it be equal to zero that corresponds to the symmetrically 
expected gusts in the y direction. Application of the expectation operator to 
(1.13) yields 


mY" (t) = mg —kY'(t), (1.14) 


where Y(t) = E[y(t)] denotes the expected trajectory. Next, equation (1.14) 
can be solved similarly to equation (1.12). 

In order to estimate the deviation of the random function y(t), first, we es- 
timate the term m(t) in (1.13) by the bounds —Co < m&(t) < Co where the 
positive constant C is determined by the chosen confidence interval estimate. 
Next, we substitute +Co in (1.13) instead of m(t) and obtain two ODE 


mY! (t) = mg — kY: (t) + Co. (1.15) 


Ultimately, we can estimate the expected trajectory and its bounds within the 
chosen confidence interval as displayed in Fig.1.2. 


6 


FIGURE 1.2: The expected trajectory (solid line) and its bounds (dashed line). 


b) Calculate the mass of the cylinder having various circular sections ex- 
pressed by the shape function S = S(x) for x € [a,b] (see Fig.1.3). It is 
assumed that the cylinder consists of the homogeneous material of the den- 
sity p. Let us divide the cylinder into small cylinders with the heights Ax; 
(i =1,2,:-- ,n). The mass of the i-th cylinder is equal to Am; = pT 8? (£;i) Azi, 
where the point €; is arbitrarily chosen in (xi, £i+1). If the latter interval is 
sufficiently small, its mass Am, is calculated quite accurately. The total mass 
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A S(x) 
| i = 
O a b x 
FIGURE 1.3: Radius S(x) of the circular section at the point x. 
of the cylinder is approximately equal to 
M, = 5 Am; = pr S? (E) Azi. (1.16) 
i=1 i=1 
Thus, the mass of the cylinder is calculated by the discrete model. 
A S(x) 
a e ' A 
O a b X 


FIGURE 1.4: Partition of the cylinder. The total mass of the cylinder is approx- 
imated by the sum of the masses of small cylinders. 


Let the value 6 = max; Az; tend to zero, i.e., [a,b] is divided into small 
intervals whose lengths become smaller and smaller. Then the sum (1.16) 
tends to a value M, which can be considered as the exact total mass of the 
cylinder. The sum Mn is called the Riemann sum in calculus. The limit M is 
called the definite integral of the function pr S?(x) on the interval (a, b): 


M= f pr’ (x) dz. (1.17) 


This integral presents the continuous model to calculate the mass. 
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Application of these two models yields a method to calculate a value (mass) 
of the complicated objects via small simple ones. 


Principle of transition: continuous © discrete. To divide a continu- 
ous object into small parts, to apply a simple formula to every part, to calculate 
the sum for all the parts and get back to the continuous object through the limit 
operation. 


One can divide an object under the study into parts in space and in time. 
Partition by the mass m = px, where p denotes the linear density, is equivalent 
to partition in space on 2. 

In mathematics, the principle of transition is the benchmark of the Rie- 
mann integral theory including multidimensional, curvilinear and surface in- 
tegrals. In continuum mechanics this principle is used in the framework of 
calculus and theory of partial differential equations. Real medium is discrete 
on the molecular level. Hence, the real mass of the cylinders is rather given by 
formula (1.16). However, summation over all molecules is practically impos- 
sible and integration is more simple than summation in this case. This does 
not concern nanotechnology when few molecules are taken in consideration. 


Remark 1.2. The above example can help to delineate the difference between 
pure mathematical and practical approaches. Assume that the radius of the 
body sections presented in Fig.1.3 is described by the Dirichlet function 


1, if x is irrational, 
D(x) = (1.18) 
0, if x is rational. 


There are rational and irrational numbers in every interval of the real axis. If 
all the points €; are irrational in the sum (1.16), then D(€;) = 1 and the total 
mass of the cylinder is equal to p X`; Ax; = p(b— a). On the other hand, if 
the points £; are rational, then D(€;) = 0. Thus, the total mass of the cylinder 
vanishes. This means that the cylinder is so porous that it does not posses a 
mass; it is between 0 and p(b — a). In this book, we take the side of applied 
mathematicians. Therefore, we do not consider such strange objects without 
any mass. However, such examples are very interesting in pure mathematics 
since they yield new mathematical topics for investigations!'. 

c) Linear models are frequently used in practice since they are simple and 
properly describe real phenomena. The definition of the linear model can be 
given in terms of the linear operator (function): 


Definition 1.1. Operator L: X — Y is called linear if for any z, £1, £2 E€ X 
and an arbitrary real (complex) number A 

i) L(x, + x2) = L(21) + L(x); 

ii) L(Ax) = AL(a). 


11Fyequently, an independently developed pure abstract theory has serious applications 
many years later after its creation. 
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Example 1.2. The function y = az is a linear operator L: R > R. 


Example 1.3. The differentiation operator 4 : f + f’ is a linear operator. 
na d d d 
This is true due to the rules 4 (f1 + fo) = “2 + ® and 4(Af) =A. 


Example 1.4. If you conduct a uniaxial tensile test on an elastic material you 
can observe that the force o per surface of unit area acting on the material is 
proportional to its deformation € = ag where AZ denote displacement of the 
specimen and £ its length. This is Hooke’s law which can be written in the 
form o = ke where the constant k depends on the material. Linear Hooke’s law 
is valid for small deformation e. Plastic transformations can occur for a large 
displacement when the coefficient k depends on e. This yields the nonlinear 
law o = k(e)e under simple loading. It is interesting to note that rubber and 
a metal are elastic materials but with different values of k. 


Cybernetics is a multidisciplinary science based on the black box when one 
considers an object which transforms inputs x € X to outputs y € Y and is 
not interested in the interior of the transformation box (Fig.1.5). 


output 


-i-— 


FIGURE 1.5: Illustration of the black box. A set of inputs and outputs determines 
the object without any study of the interior structure. 


input 


The rule of the transformation K : x — y is interesting in cybernetics 
(roughly speaking, software structure of K, not its hardware). That is why 
this object is called the black box. Such an approach is frequently used in the 
theory of signals. A signal x(t) depending on time t is transformed into y(t) 
by the integral operator: 


T 
(Ka) (t) := y(t) = | k(7, t)a(r) dr. (1.19) 


One can check that the operator (1.19) is linear. Let the kernel have the special 
structure k(7,t) = k(t — t). Then, the operator 


T 
(Kx)(t) := y(t) = f k(T — t)a(r) dr (1.20) 


is called the convolution of k(T) and x(r). The transformation K : x —> y can 
be non-linear [23] 


T 
(K2)(t) := y(t) = f k(T,t,£(T)) dr, (1.21) 


where k(r, t,x) is a function of three variables. 
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1.3 Stability of models 


One of the most important characteristic of the model is its stability. Me- 
chanical illustration of stability is displayed in Fig.1.6. Stability can be dis- 
cussed on a parameter and in time. 


a) b) 


FIGURE 1.6: Stable (a) and not stable (b) locations of a ball 


Consider the function (1.5) which is continuously dependent on the initial 
velocity v (parameter of the model) in the closed interval [0,7]. This means 
that for any t € [0,7] and for sufficiently small changes of vo the function 
(1.5) has small increments. More precisely, let vı and vg be two such values 
of the initial velocity that |v, — v| < £ for any fixed positive number e. Let 
yi(t) and y2(t) be functions of the form (1.5) with the velocities vı and v2, 
respectively. Then 


ma lyi(t) — y2(t)| = nag |vit — vat] = T|v1 — v2| < Te. 
This inequality implies that for sufficiently small € > 0, the maximal value 
of the difference |yi(t) — y2(t)| for t € [0, T] is also small, i.e., the discussed 
problem is stable for vo. 

On the other hand, the time T can be a large number, it can even be equal 
to +oo. Then max; 0, +0] [yi (t) — y2(t)| = +00 for |vı — v2| < 0.000001. This 
means that the problem is unstable in v. 

A similar situation takes place for weather prediction (see the footnote 
on page 16 with the reference to Lorenz’s paper and the advanced book [6]). 
Equations of the air flow on the Earth scale can be stable. Let us say for 
simplicity that a function y(t) describing the air flow satisfies inequality 


mane [on (0) = ya] < ele. 
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Here, the constant c(T) depends on the prediction time T. The value c(T’) can 
be bounded for T which lasts for about a week, what yields well estimations 
for y(t). However, c(T) abruptly increases for greater T in such a way that 
c(T)e also increases. The stable process becomes unpredictable. 

In main, stability is treated as T tends to infinity. Linear homogeneous 
ODE with constant coefficients have solutions presented as a sum of modes. 
Each mode has the form y(t) = P(t) exp(At) where P(t) is a polynomial, A = 
a+i@ is a complex constant and i denotes the imaginary unit. A mode is called 
stable if it is bounded as t + +00. Euler’s formula exp(At) = e®(cos Bt + 
isin Bt) yields 


Theorem 1.1. Let P(t) be a non-constant polynomial. A mode y(t) = 
P(t) exp(At) is stable if and only ifReA=a< 0. 


Formally, one can divide models into stable and unstable ones. But a physi- 
cist will not permit us to do it referring to the principle of stability: 


Principle of stability. A mathematical model must be stable. 


The study of the unstable problem should not give reasonable results since 
small exterior changes or small changes in the initial data can yield large devi- 
ation in the obtained results. But at the beginning of the study no one knows 
whether the model is stable. Therefore, one should first develop a model, in- 
vestigate it and only after its verification it is possible to claim its stability. It 
follows from the above examples that the same model can be stable and unsta- 
ble. It is worth investigating the conditions under which the problem is stable. 


Concluding remarks and further reading concerning this section. 
Not everyone accepts the principle of stability in such a form and investigates 
unstable processes following the principle by Dirac stated on page 126. This 
can be reasonable, for instance, in biological and chemical pattern formations 
(22, 33]. In the theory of turbulence, the amplitude of the average velocity can 
be considered when the local velocity is unpredictable [2]. 

Various mathematical approaches to stability are discussed in courses of 
ODE and PDE (see Sec.4.3, Sec.4.8 and [5], [41], [83], [22]). The numerical 
stability is discussed in Sec.8.6.1 ([35], [52]). 


1.4 Dimension, units, and scaling 
1.4.1 Dimensional analysis 


Mathematical models contain measurable values. One cannot say that a 
dimension value is large or small because it depends on the units in which 
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it is expressed. For example, the length L can be measured in meters (m), 
centimeters (cm), miles etc. The length L will be called an abstract dimension 
and the length m a special one. One can treat L as a variable which takes 
special values m, cm, etc., which can be written as L = m. Everybody may 
introduce his/her own dimension unit. It is only necessary to relate it to 
the accepted units. The following three dimensions of units are accepted in 
physics and technology: SI, CGS and English units. Besides the length L 
we introduce the abstract units of mass M and time T. The velocity has 
the abstract unit & and the special one with L = m and T = sec where 
sec denotes second. It is convenient to introduce brackets for the units. For 
instance, [op] = ML? means that the density p is considered in dimension 
units and measured in mass per volume. Similarly, [53] = kg m~3 means 
that the density is measured in kilo per spatial meter and is equal to 5 in these 
units. The expression [p’] = 1 denotes that p’ is measured in dimensionless 
units. 

Consider the example from Sec.1.1.1 in the abstract units [m] = M, 
t] = T, [y] = L, vo = LT, [g] = LT~?. We are interested in the units 
which appear in equation (1.5). The first term g on the right-hand side has 


the dimension [2] = [g|[t}? = LT~?T? = L. The second term vot on the 
right-hand side has the dimension [vot] = [vo][t] = LT~'T = L. Hence, the 
right-hand side contains the sum of two values of the same dimension L, that 
confirms the proper dimensional usage of the expression ue + vot. The left- 
hand side of the equation also has the dimension L. Everything is alright. 
Checking the dimensions yields a simple way to check an equation. For in- 
stance, if one develops a model and gets equations L = L + M or T = LT”, 
this gives a signal that something is wrong because a meter cannot be com- 
pared to a kilo. 


Principle of units. It is forbidden to add camels to tractors. But it is 
possible to multiply them. 


The units of derivatives can be obtained from the definition of the deriva- 
tive through the limit of the difference quotients. For instance, [y] = L implies 
[y'] = LT} and [y”] = LT~?. We now proceed with considering the units of 
equation (1.12) 


[my"] = [mg] + [ko] = [m][y"] = [m][9] + [Alle] 
& MLT? = MLT? + [k] LTt. 


The principle of units holds if and only if [k] LT} = MLT~?. This implies 
that the unit of the air resistance, the coefficient k, must be equal to MTŁ. 
The coefficient k can be estimated in experiments. Let M = kg and T = sec. 
In order to determine k one has to present the experimental results in such a 
way that k will be measured in kilo per second. 

Dimensional analysis is applied not only to check a dimensional equation 
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but can suggest a type of such an equation. The main fundamental theorem of 
dimensional analysis is the Pi theorem [41]. We formulate it in the following 
reduced form 


Theorem 1.2. Let an equation (E) relate N dimensional quantities. There 
exists such an equation (e) which contains n dimensionless quantities and 

i) equations (e) and (E) are equivalent; 

i)n <N. 


The full form of the Pi theorem is based on the dimension matrix. We 
explain this notation on the equation (1.5) written in the form 


g? 
F(y,t,9,%0) =y- “> — vot = 0. (1.22) 
The dimensions are given by [y] = L, [t] = T, [vo] = LT} and [g] = LT~?. 
The dimension matrix A is defined by the powers occurred at these equations. 
It is convenient to write A in the form of the table 


y t | vw | g 
L I 0 I I 
T 0 1] —1 | —2 


Two dimensionless values can be introduced 


ne ee nee (1.23) 
y 
The value g and vo are expressed through other values by equations 


oe: sY 
gaT Vona (1.24) 
Substitution of (1.24) into (1.22) yields 


1— om = 0. (1.25) 


An experimental verification of (1.22) can be done by measurement of y, t, vo 
and g in arbitrary units. Further, mı and 7 are calculated with (1.23). At the 
end, the value 1 — 4 — m2 has to be compared with zero, i.e., equation (1.25) 
has to be checked. 


1.4.2 Scaling 


Each dimensional parameter has its domain of values. Estimations of this 
domain are useful in the introduction of the dimensionless units by scaling. 
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The parameters of equation (1.12) can be made dimensionless in the following 
way 


y = 9, t= TË, m= mon, (1.26) 


where £, T, mo are characteristic values of the length, time and mass, respec- 
tively. For instance, if a book falls from a table, one can assume £ = lcm, 
T = lsek, mo = 0.1kg. If a parachutist falls from an plane, one can take 
£ = 700m, T = Imin, mo = 10kg. The choice of the scale can be arbitrary. 
However, it is better to take reasonable values in order to work with “usual” 
numbers, to avoid divisions into small numbers following the principle of the 
balanced computations. 
The derivatives y’ and y” after the substitutions (1.26) are transformed as 
follows 7 j zs 
a dy _ fay (1.27) 
dt rd dt? 71? dt? 
Then equation (1.12) becomes 


i’ + kil = 5, (1.28) 


where k = kr denotes a dimensionless coefficient. This can be proved in the 
following way 7 
[k] = [k] TM} = MTTM =1. 


Concluding remarks and further reading concerning this section. 
A short introduction to the unit theory and scaling is presented above. Further 
extensions and applications can be found in [41]. 


Exercises 


1. Describe the trajectory of the falling object thrown in a horizontal di- 
rection. 


2. Describe the trajectory of the falling object thrown in a horizontal direc- 
tion when the land is determined by a surface, for instance by a linear 
function y = a1£1ı + a2£2 + b. 


3. Fig.1.1 presents a film frame of the falling object. Prepare an animation 
of the whole process. Calculate online the velocity and the acceleration 
of the point. Take into account the air resistance. 


Useful Mathematica operators: Animate, Manipulate, Slider, 
TrackedSymbols, DynamicModule. For creating animation in MAT- 
LAB refer to Example Box 8.1 on page 175. 
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As an advanced exercise, you can create slider manipulator in MATLAB. 
Refer to the documentation and investigate the uicontrol operator. 


. Calculate numerically the sums XZ; +; for m = 1,2,3,4,.... 


n=l n™ 
Think about the outputs Zeta[3] and Zeta[3.] in Mathematica for 
m = 3. 


. Let hot coffees and cold cream in Example 1.1 be replaced by cold water 


with ice and juice at room temperature, respectively. Who drinks colder 
water? 


. Calculate the mass of the circular cone of height h and radius r with 


constant density p. 


. Are the following operators linear? 


A linear function y = ax + b from R on R; 
function Y = AX from R” in R”, where A is a n x n matrix; 


integral operator [>~ k(r, t)x(t) dr. 


. Present a graphical interpretation of the stability of the process (1.5) in 


time when the initial velocity vo is measured with an error of 5%. 


